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ABSTRACT 

These  lecture  notes  present  the  basic  theory  and  methods  for  the  Direct  Simulation  Monte  Carlo  (DSMC) 
algorithm.  Some  of  the  open  challenges  in  the  treatment  of  complex,  multi-scale  flows  are  also  discussed. 


1.0  INTRODUCTION 
1.1  Dynamics  of  Dilute  Gases 

Ligure  1  illustrates  the  large  range  of  physical  scales  for  dilute  gases,  ranging  from  the  quantum  scale  to 
the  hydrodynamic  scale.  At  the  smallest  scales,  the  details  of  the  intermolecular  interaction  are  important; 
the  corresponding  time  scale  is  the  duration  time  of  a  collision,  which  is  on  the  order  of  femtoseconds.  At 
the  opposite  end,  when  the  gas  is  well-represented  by  continuum  hydrodynamic  fields  (density, 
temperature,  etc.)  then  a  description  based  on  partial  differential  equations,  such  as  the  Navier-Stokes 
equations,  is  accurate. 
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At  the  intermediate,  kinetic  scale  the  discrete  molecular  nature  of  the  gas  has  a  significant  effect  on  the 
dynamics.  The  relevant  length  and  time  scales  in  this  regime  are  the  mean  free  path  and  mean  free  time 
and  they  establish  all  the  transport  properties.  At  this  scale  the  continuum  approximations  for  transport, 
such  as  the  Lourier  law  for  heat  flux,  are  not  accurate  so  standard  Computational  Lluid  Dynamics  (CLD) 
approaches  are  not  suitable.  On  the  other  hand,  the  collisions  themselves  may  accurately  be  treated  using  a 
simplified  classical  approximation  or  as  stochastic  events  with  specified  rates. 
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1.2  Direct  Simulation  Monte  Carlo 

Direct  Simulation  Monte  Carlo  (DSMC)  is  a  molecular  scheme  designed  for  the  simulation  of  gases  at  the 
kinetic  scale.  [1]  DSMC  was  developed  by  Graeme  Bird  in  the  late  1960’s  for  kinetic  scale  simulations 
and  was  quickly  adopted  by  the  aerospace  engineering  community  because  the  method  is  accurate  for 
flows  with  high  Knudsen  number  (Kn),  the  ratio  of  mean  free  path  to  system  length.  For  example,  during 
the  initial  re-entry  flight  in  the  upper  atmosphere  the  mean  free  path  in  the  air  is  comparable  to  the  size  of 
the  vehicle.  Although  DSMC  was  originally  motivated  by  these  “Space  Age”  applications,  the  method  is 
now  used  in  a  wide  range  of  problems  in  physics,  chemistry,  and  engineering.  Examples  range  from 
micron-scale  flows  (which  are  also  high  Knudsen  number)  to  granular  gases  to  lunar  atmospheres.  The 
algorithm  has  evolved  as  well,  with  improvements  to  the  numerical  accuracy  and  efficiency  as  well  as 
extensions  to  complex  chemistry  and  even  to  dense  gases.  Nearly  50  years  after  its  introduction  DSMC 
remains  the  dominant  numerical  method  for  molecular  simulations  of  dilute  gases. 

1.3  Molecular  Dynamics  and  DSMC 

Molecular  Dynamics  (MD)  is  another  popular  numerical  approach  for  molecular  simulations  but  it  is 
normally  applied  to  liquids.  [2]  In  MD  the  exact  molecular  motion  is  calculated  from  the  inter-molecular 
forces  among  the  particles.  Although  this  may  be  done  efficiently  for  liquids,  in  which  the  molecules 
directly  interact  with  only  a  few  neighbouring  molecules,  in  a  dilute  gas  a  molecule  has  thousands  of 
potential  collision  partners.  Direct  Simulation  Monte  Carlo  (DSMC)  overcomes  the  inefficiency  of  MD  by 
replacing  the  deterministic  motion  with  a  stochastic  approximation  for  the  collision  process.  DSMC  is  able 
to  advance  in  time  steps  comparable  to  the  mean  free  time  between  collisions  yet  remain  accurate  at  the 
level  of  the  Boltzmann  equation.  Unlike  MD,  DSMC  is  always  numerically  stable  regardless  of  time  step. 


2.0  DSMC  METHOD 

2.1  Outline  of  DSMC 

Before  describing  the  details  of  how  DSMC  works,  let  us  quickly  outline  the  algorithm’s  framework. 
DSMC  is  a  particle-based  scheme  so  a  typical  calculation  initializes  the  desired  geometry  with  boundary 
conditions  and  fills  the  computational  volume  with  random  particles.  At  each  time  step  all  particles  move 
ballistically  according  to  their  assigned  velocity;  in  DSMC  collisions  are  independent  of  these  trajectories. 
Any  particles  reaching  a  boundary  are  processed  according  to  the  imposed  conditions,  such  as  randomly 
assigning  a  new  velocity  to  a  particle  that  strikes  a  thermal  wall.  If  there  are  open  boundaries  (e.g.,  wind- 
tunnel  configuration)  then  particles  are  generated  as  inflow  and  removed  if  they  exit  through  the  boundary. 
Quantities  of  interest  in  a  simulation  are  either  flow  properties  (e.g.,  temperature  field  in  the  flow)  or 
surface  quantities  (e.g.,  drag  and  lift  for  a  vehicle)  and  these  are  measured  by  statistical  samples  of 
particles  contained  in  a  volume  or  striking  a  surface.  Finally,  DSMC  randomly  selects  and  evaluates 
collisions  among  the  particles.  This  last  step  is  the  heart  of  the  algorithm  so  it  will  be  described  in  some 
detail. 

2.2  Random  Numbers 

As  its  name  implies,  Direct  Simulation  Monte  Carlo  uses  random  numbers.  Unlike  other  Monte  Carlo 
schemes,  such  as  Metropolis  MC  or  Quantum  MC,  DSMC  uses  a  wide  variety  of  probability  distributions 
for  different  purposes.  For  example,  to  initialize  particles  in  a  volume  we  might  first  determine  the  number 
of  particles  to  insert  by  choosing  the  value  from  a  Poisson  distribution.  If  the  particles  are  distributed 
uniformly  then  their  coordinates  are  generated  from  the  uniform  distribution;  their  velocity  components 
are  usually  generated  from  the  Maxwell-Boltzmann  distribution,  which  is  equivalent  to  a  normal  or 
Guassian  distribution.  As  described  below,  at  various  stages  of  the  advection  and  collision  processes  more 
random  variables  are  used  from  various  distributions. 
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Given  the  importance  of  the  Monte  Carlo  elements  in  DSMC  it  is  essential  to  use  a  high-quality  random 
number  generator.  Simple  generators  that  were  adequate  a  decade  ago  are  obsolete  today  because,  with 
increased  computer  power,  modem  DSMC  calculations  can  cycle  through  the  entire  period  of  a  simple 
generator  in  a  matter  of  seconds.  Although  DSMC  generates  random  values  from  a  variety  of  distributions, 
the  program  only  needs  one  generator  that  produces  the  uniform  distribution  in  (0,1);  all  other 
distributions  are  obtained  by  transformation  (usually  inversion  of  the  cumulative  distribution)  or  iteration 
(typically  an  acceptance  /  rejection  scheme).  [3]  There  are  several  modem  generators  for  the  uniform 
distribution  to  choose  from;  the  author’s  current  favourite  is  the  Mersenne  Twister.  [4]  One  final  note 
regarding  the  use  of  random  number  generators:  for  the  purpose  of  debugging  the  initial  state  (or  seed) 
should  be  recorded  so  that  a  calculation  that  fails  can  be  repeated  exactly. 

2.3  Boundary  Conditions 

Once  the  particles  have  been  initialized  the  DSMC  calculation  advances  in  time  steps  At,  alternating 
between  ballistic  motion  of  the  particles  and  collisions  among  the  particles.  Typically  the  particles  simply 
are  displaced  by  a  distance  Ar  =  v  At  but  if  there  is  a  body  force,  such  as  gravity,  then  the  motion  is 
slightly  more  complicated.  The  particles  move  without  interaction  and  can  even  overlap;  the  only  place  in 
DSMC  where  the  particles’  cross-section  is  used  is  in  determining  the  collision  rate.  For  transient  flows, 
on  the  first  time  step  one  should  use  V2  At  (Strang  splitting)  to  maintain  temporal  accuracy.  If  measuring 
non-conserved  variables  (e.g.,  fluxes)  then  one  should  also  time-center  the  sampling  (half  move,  sample, 
half  move,  collisions)  for  all  steps. 

During  their  motion  some  of  the  particles  will  reach  boundaries,  either  at  the  periphery  of  the 
computational  domain  or  at  the  surfaces  of  objects  within  the  domain.  One  of  the  strengths  of  DSMC  is 
the  ease  with  which  boundary  conditions  may  be  imposed.  For  periodic  boundaries  a  particle  crosses  one 
boundary  and  re-enters  at  the  periodic  mirror  boundary.  For  specular  boundaries  a  particle  reflects  off  the 
surface  with  the  perpendicular  component  of  velocity  reversed.  At  a  thermal  wall  a  particle’s  velocity 
components  are  reset  according  to  the  biased-Maxwellian  distribution  for  the  normal  component 


v^e-m.v\/2kTm 


and  the  standard  Maxwell-Boltzmann  distribution  for  the  parallel  components, 


Both  of  these  distributions  are  easily  generated  by  inversion,  the  former  from  an  exponential  distribution 
and  the  latter  from  a  normal  distribution.  A  surface  may  not  be  fully  accommodating,  for  example  in  the 
Maxwell  model  a  random  fraction,  a,  are  thermalized  while  the  rest  reflect  specularly. 

Inflow/outflow  boundary  conditions  are  commonly  treated  as  a  reservoir  with  given  density,  fluid  velocity, 
temperature.  Particles  in  the  main  system  are  removed  if  they  cross  the  boundary  into  the  reservoir. 
Particles  are  injected  from  reservoir  to  main  system  by  either:  Surface  generator  -  From  the  number  flux 
determine  number  to  be  injected  during  a  time  step  and  then  generate  particle  velocities  from  surface 
distribution  (e.g.,  inflow  Maxwellian);  Volume  generator  -  Fill  a  “ghost  cell”  with  particles  before  the 
ballistic  move  and  discard  any  that  do  not  cross  the  boundary  into  the  main  system  during  the  move  phase. 
Efficient  schemes  are  available  for  both  types  of  generators.  [5,6] 
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ORGANIZATION 


2.4  Collision  Selection 

To  evaluate  collisions  the  domain  is  partitioned  into  cells  and  the  particles  are  sorted  into  them;  the 
algorithm  then  loops  over  cells  evaluating  collisions  independently  in  each  cell.  A  maximum  collision 
rate,  7?max  =  R (N,  V ,  o,  max{vr}),  is  calculated  in  each  cell  based  on  the  number  of  particles,  the  cell 
volume,  the  collisional  cross-section  for  the  particles,  and  the  (estimated)  maximum  relative  speed.  The 
explicit  form  for  this  function  depends  on  the  interaction  potential;  for  example,  for  hard  spheres, 


Ni.N  —  iters- 

2V 


where  a  =  nd2.  The  mean  collision  rate  is 


m — 

since  N  is  Poisson  distributed.  The  number  of  attempted  collisions  during  a  time  step  is  7?max  At ;  for  each 
attempted  collision  a  random  pair  of  particles  is  selected  within  the  cell.  Given  the  actual  relative  velocity 
for  the  selected  pair,  they  are  accepted  as  collision  partners  with  probability  R(N,  V ,  a,  vr)/i?max.  This  not 
only  accepts  pairs  with  the  correct  probability  (pairs  with  higher  relative  velocity  are  more  likely  to 
collide)  but  also  accepts  them  with  the  correct  mean  collision  rate. 

2.5  Collision  Evaluation 

Once  a  pair  is  accepted  for  a  collision  all  that  remains  is  the  calculation  of  the  particles’  post-collision 
velocities  (if  there  are  also  internal  degrees  of  freedom  or  chemistry  then  there  are  addition  steps).  By 
conservation  of  momentum  the  center-of-mass  velocity  is  constant;  for  particles  of  equal  mass, 
conservation  of  energy  implies  that  the  magnitude  of  the  relative  velocity  is  unchanged.  These  constraints 
fix  four  of  the  six  unknowns;  the  remaining  two  are  set  by  the  direction  angle  for  the  relative  velocity. 
Because  DSMC  does  not  use  the  actual  particle  trajectories  to  evaluate  collisions  this  direction  angle  is 
selected  at  random.  For  hard-sphere  collisions  the  direction  angle  is  chosen  as  uniformly  distributed  in  the 
unit  sphere  by  inversion.  [3] 

2.6  DSMC  “Parliament” 

In  DSMC  the  number  of  simulation  particles  (“simulators”)  is  typically  a  small  fraction  of  the  number  of 
physical  molecules  with  each  simulator  representing  Aef  physical  molecules.  DSMC’s  dynamics  is  correct 
if:  a)  The  DSMC  simulators  are  an  unbiased  sample  of  the  physical  population  (unbiased  parliament);  b) 
Collision  rate  is  increased  by  Ne f  so  the  number  of  collisions  per  unit  time  for  a  simulator  is  same  as  for  a 
physical  molecule;  c)  In  sampling,  each  simulator  counts  as  Ne f  physical  molecules.  Early  DSMC 
implementations  used  a  different  (unpopular)  representation,  rescaling  the  simulator  diameter  and  mass  to 
maintain  the  same  physical  mean  free  path  and  mass  density.  While  increasing  NG f  and  lowering  the 
number  of  simulators  speeds  the  calculation  it  also  decreases  the  accuracy.  Empirically  it  has  been  found 
that  the  accuracy  of  DSMC  goes  as  l/N;  for  traditional  DSMC  keeping  about  N=  20  particles  per  collision 
cell  is  the  rule-of-thumb.  [7] 

2.7  Ballistic  and  Collisional  Transport 

By  their  ballistic  motion  particles  carry  mass,  momentum  and  energy.  In  a  dilute  gas,  this  is  the  only 
physical  mechanism  for  transport.  Yet  in  DSMC,  momentum  and  energy  are  also  transported  by  the 
collisions.  The  larger  the  collision  cell,  the  more  collisional  transport  because  of  the  greater  average 
separation  between  particle  pairs.  The  collisional  transport  may  be  calculated  by  Green-Kubo  theory, 
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which  finds  that  this  error  is  quadratic  in  cell  size  and  time  step.  [8,9]  To  minimize  this  collisional 
transport  the  rule-of-thumb  for  traditional  DSMC  is  to  limit  the  cell  size  to  a  fraction  of  a  mean  free  path 
and  the  time  step  to  a  fraction  of  a  mean  collision  time.  New  DSMC  implementations  minimize  collisional 
transport  error  by  choosing  the  collision  partner  as  the  nearest  particle  in  the  cell.  To  avoid  bias  due  to  re¬ 
collisions,  a  particle  pair  is  not  allowed  to  collide  twice  (choose  next-nearest  particle  instead).  Two 
common  implementations  are:  Transient  Adaptive  Sub-cells  (introduced  by  Bird);  and  Virtual  Sub-cells 
(introduced  by  LeBeau,  et  al).  [10] 


3.0  HYDRODYNAMIC  FLUCTUATIONS  AND  STATISTICAL  ERROR 
3.1  Fluctuations  in  DSMC 

Even  at  thermodynamic  equilibrium,  the  state  of  the  system  in  DSMC  has  random  variations  as  particles 
move  and  collide.  It  is  important  to  note  that  these  hydrodynamic  fluctuations  in  density,  temperature,  etc. 
have  nothing  to  do  with  the  Monte  Carlo  aspect  of  DSMC.  Molecular  Dynamics  is  a  deterministic 
algorithm  and  it  has  precisely  the  same  fluctuations.  The  random  variations  in  particle-based  schemes  arise 
from  the  complex  ergotic  mixing  dynamics  due  to  the  non-linear  interactions. 

The  variance  of  fluctuations  in  DSMC  is  exact  at  equilibrium;  particles  are  uniformly  distributed  in 
position  and  their  velocities  are  Maxwell-Boltzmann  distributed.  The  spatial  and  temporal  correlations  of 
fluctuations  in  DSMC  are  also  correct  at  hydrodynamic  scales,  as  demonstrated  by  the  agreement  found 
with  Landau-Lifshitz  fluctuating  hydrodynamics,  and  beyond.  [11]  Surprisingly,  this  agreement  extends  to 
strongly  non-equilibrium  scenarios  such  as  severe  temperature  gradients  and  strong  shearing  flows.  [12] 

If  the  phenomenon  of  interest  has  dynamics  that  depend  on  the  accurate  representation  of  fluctuations, 
such  as  Brownian  motion,  then  the  fluctuations  in  DSMC  are  a  positive  feature.  However,  for  most 
engineering  applications  these  fluctuations  are  an  annoyance  because  of  the  amount  of  sampling  that  has 
to  be  done  in  order  to  extract  accurate  answers.  Furthermore,  when  each  DSMC  simulation  particle 
represents  a  large  number  of  physical  molecules,  as  is  typical  in  engineering  calculations,  the  variance  of 
fluctuations  is  magnified  by  this  ratio.  For  this  reason  fluctuations  are  unphysically  large  when  Nef>  1. 


3.2  Error  Bars  in  DSMC 

Fortunately,  because  the  variances  of  fluctuations  in  DSMC  are  known  from  statistical  mechanics,  we  may 
estimate  a  priori  the  amount  of  sampling  needed  to  obtain  the  desired  accuracy.  Furthermore,  since  the 
non-equilibrium  corrections  are  small,  even  when  gradients  are  large,  we  may  use  the  simpler  results  for 
equilibrium  fluctuations  for  estimating  error  bars.  [13]  For  example,  the  fractional  error  in  the  sample 
mean  of  the  x-component  of  fluid  velocity  is 


E 


u 


_  i _ L 

k  I  ~  y[SN  Ma 


where  S  is  the  number  of  samples,  N  is  the  average  number  of  DSMC  particles  in  the  sampling  volume, 
and  Ma  is  the  Mach  number  of  the  flow.  We  immediately  see  that  in  aerospace  applications,  where  the 
Mach  number  is  large,  only  a  small  number  of  samples  is  required.  In  microscopic  flows,  where  Ma  «  1, 
the  fractional  error  is  much  greater  and  far  more  sampling  is  required,  especially  since  the  error  bar 
magnitude  decreases  as  S~m.  Similar  results  are  found  for  other  hydrodynamic  quantities  and  these  results 
apply  equally  to  other  particle-based  schemes,  such  as  Molecular  Dynamics. 
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3.3  Variance  Reduction 

Attempts  to  reduce  the  variance  in  DSMC  have  had  mixed  success.  Some  approaches  are  effective  but  at 
the  expense  of  introducing  significant  complexity  to  the  calculations  of  collisions,  advection,  boundary 
conditions,  etc.  Other  approaches  are  simple  to  implement  but  are  not  always  numerically  stable.  Worst  of 
all  are  those  approaches  that  inadvertently  introduce  statistical  errors  that  bias  the  sampling.  For  example, 
if  the  cells  into  which  particles  are  partitioned  are  dynamically  resized  (for  example,  to  use  larger  cell 
volumes  in  regions  with  few  simulation  particles)  then  the  cell  volume  is  a  random  variable,  possibly 
correlated  to  other  variables.  In  general  <1/V>  f  H<V>  so  calculations  involving  density  (e.g.,  collision 
rate)  will  have  a  statistical  bias  depending  on  the  probability  distribution  for  volume.  Nevertheless,  there 
has  recently  been  promising  progress  in  variance  reduction  for  DSMC  (and  its  variants)  and  this  remains 
an  active  and  important  topic  for  research.  [14,  15] 

3.1  Thermal  Fluctuations  and  Brownian  Dynamics 

Returning  to  the  positive  viewpoint  of  fluctuations,  there  are  a  number  of  important  applications  where  the 
stochastic  nature  of  a  fluid  at  microscopic  scales  is  central  to  the  phenomenon  of  interest.  The  accurate 
modelling  of  multi-scale  fluid  phenomena  often  requires  the  correct  representation  of  the  spatial  and 
temporal  spectra  of  fluctuations,  specifically  when  studying  systems  where  the  microscopic  stochastics 
drive  macroscopic  phenomena.  Some  examples  in  which  spontaneous  fluctuations  play  an  important  role 
in  hydrodynamics  include  the  breakup  of  droplets  in  jets,  Rayleigh-Bemard  convection,  Kolmogorov  flow, 
and  Rayleigh-Taylor  mixing.  Thermal  fluctuations  are  also  important  in  many  chemical  phenomena,  such 
as  combustion  and  explosive  detonation  as  well  as  the  propagation  of  reaction  fronts.  Finally,  micro¬ 
machines  have  been  proposed  that  make  use  of  spontaneous  fluctuations  as  an  essential  part  of  the  design. 
In  fact  many  biological  mechanism  operating  at  the  cellular  level  employ  such  “Brownian  motors.”  In 
Directing  Matter  and  Energy:  Five  Challenges  for  Science  and  the  Imagination ,  a  report  by  the  Basic 
Energy  Sciences  Advisory  Committee,  the  authors  state  that,  "evolution  has  embraced  stochastic 
fluctuations  and  often  relies  on  them  for  the  functionality  of  the  system.  This  suggests  an  interesting 
design  principle  that  humans  have  not  yet  learned  to  use.  Exploitation  of  statistical  fluctuations  may  well 
be  essential  to  accomplish  some  of  the  more  exotic  tasks  living  systems  are  able  to  perform... Realizing  the 
promise  of  nanoscience  requires  that  we  deal  with  non-equilibrium  and  fluctuations."  Because  DSMC 
simulations  include  spontaneous  fluctuations  with  the  correct  physical  spectrum  (when  Ne f  =  1)  the  method 
is  ideally  suited  for  these  types  of  problems  in  dilute  gases. 


4.0  ALGORITHM  REFINEMENT 
4.1  Multi-Scale  Modeling 

Algorithm  Refinement  (AR)  is  an  important  new  approach  for  modeling  multiscale  fluid  phenomena. 
Based  on  the  framework  of  mesh  refinement,  AR  combines  two  or  more  algorithms  that  use  different 
descriptions  of  the  flow  at  different  scales.  Such  hybrid  schemes  typically  couple  structurally  different 
computational  schemes  such  as  particle-based  molecular  simulations  with  continuum  partial  differential 
equation  (PDE)  solvers.  The  general  idea  is  to  perform  detailed  calculations  using  an  accurate  but 
expensive  algorithm  in  a  small  region  (or  for  a  short  time),  and  couple  this  computation  to  a  simpler,  less 
expensive  method  applied  elsewhere. 

Algorithm  Refinement  addresses  several  difficulties  associated  with  multi-scale  numerical  modelling. 
First,  the  wide  range  of  salient  scales  makes  it  impractical  or  impossible,  even  with  the  next  generation  of 
supercomputers,  to  perform  full  system  calculations  using  algorithms  that  capture  the  correct  physics  at 
the  finest  resolution.  Second,  numerical  approaches  that  rely  on  continuum  models  and  constitutive 
closure  relations  are  limited  by  the  underlying  assumptions  (which  are  often  phenomenological  and 
sometimes  not  well  understood),  leading  to  erroneous  predictions  if  used  in  regimes  where  these 
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approximations  break  down.  Finally,  mathematical  models  use  disparately  different  representations  for 
molecular,  mesoscopic,  or  macroscopic  scales,  and  the  corresponding  algorithms  echo  this  disparity. 
Although  the  mathematical  models  and  corresponding  numerical  methodologies  used  in  hybrid  schemes 
are  often  mature  techniques,  problems  can  arise  because  the  various  methods  are  usually  developed 
independently  with  little  consideration  to  their  coupling  to  other  models  or  methods  appropriate  for  other 
scales.  Naive  attempts  to  reduce  the  computational  expense  of  molecular  simulations  by  combining  them 
with  simple  continuum  models  have  often  resulted  in  ad-hoc  constructions  that  are  numerically  inaccurate, 
unstable,  and  inefficient  because  the  two  algorithms  were  simply  cobbled  together. 

4.2  Algorithm  Refinement  with  DSMC 

The  coupling  of  DSMC  with  a  continuum  CFD  solver  was  introduced  over  20  years  ago  [16]  and  many 
implementations  have  been  presented  in  the  literature.  Most  of  these  hybrids  have  been  aerospace 
engineering  codes;  as  described  in  the  previous  section,  spontaneous  fluctuations  in  DSMC  are  an 
annoyance  for  those  applications.  On  the  other  hand,  the  group  at  Lawrence  Berkeley  lab  has  been 
developing  DSMC/CFD  hybrids  with  an  eye  towards  complex,  stochastic  flows.  Figure  2  schematically 
illustrates  the  formulation  of  an  AR  hybrid  for  the  Brownian  motion  of  a  sphere. 


Brownian  particle 


Molecular  simulation 
of  solvent  fluid 


Interface 


Continuum  simulation 
of  solvent  fluid 


Figure  2:  Algorithm  Refinement  hybrid  coupling  a  DSMC  calculation  near  the 
Brownian  particle  and  a  continuum  PDE  calculation  in  the  bulk  of  the  fluid. 


Our  research  has  been  directed  along  three  avenues:  First,  our  original  AR  scheme  [17]  was  limited  to 
dilute  gases  using  DSMC  but  we  have  developed  more  advanced  stochastic  particle  schemes  for  non-ideal 
fluids.  These  are  discussed  in  Section  5.  Second,  perfecting  the  coupling  of  particle  and  PDE  schemes  is 
challenging  due  to  fluctuations.  Though  our  coupling  is  still  not  entirely  seamless,  we  have  a  better 
understanding  as  to  the  physical  origin  of  the  mismatch.  Finally,  in  the  past  we  used  a  deterministic, 
explicit  scheme  for  the  full  Navier-Stokes  equations  but  we  now  have  efficient,  accurate  stochastic  PDE 
schemes  for  the  Landau-Lifshitz  fluctuating  hydrodynamic  equations. [18, 19] 

4.3  Hybrids  with  Stochastic  PDE  Solvers 

While  a  stochastic  PDE  solver  is  not  difficult  to  implement,  the  question  naturally  arises:  Is  it  necessary  to 
include  accurate  hydrodynamic  fluctuations  in  the  continuum  region  given  that  the  particle  region,  which 
surrounds  the  Brownian  particle,  already  has  fluctuations?  The  answer  is  a  resounding  yes,  as  we 
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discovered  in  two  basic  examples.  [20]  First,  the  velocity  auto-correlation  function  for  the  Brownian 
particle  is  accurately  reproduced  by  an  AR  hybrid  with  stochastic  PDEs  but  not  with  deterministic  PDEs. 
The  second  example  is  the  “adiabatic  piston”,  illustrated  in  Figure  3.  Molecules  collide  specularly  with  the 
massive  piston,  which  then  performs  an  asymmetric  Brownian  walk  until  the  gases  on  the  two  sides  reach 
mechanical  and  thermal  equilibrium  (i.e.,  equal  pressure  and  temperature).  In  Figure  4  we  clearly  see  that 
the  equilibration  of  the  piston  position  with  time  is  incorrect  for  the  deterministic  hybrid,  even  though  the 
particle  region  is  centred  adaptively  to  follow  the  piston  position. 


Figure  3:  Adiabatic  Piston  problem.  Schematic  (left);  AR  hybrid  snap-shot  (right). 


Figure  4:  Time-relaxation  of  the  position  of  the  adiabatic  piston  in  AR  hybrid  simulations. 
Initial  mechanical  equilibrium  (left);  initial  mechanical  non-equilibrium  (right). 


5.0  NON-IDEAL  GASES 
5.1  Equation  of  State  in  DSMC 

In  DSMC  the  various  collision  models  allow  for  realistic  representations  of  internal  degrees  of  freedom, 
which  yield  accurate  heat  capacities,  and  of  transport  properties,  such  as  the  viscosity  and  thermal 
conductivity.  However,  the  equation  of  state  in  standard  DSMC  is  always  that  of  an  ideal  gas,  P  =  Nk  T/V. 
Fundamentally,  pressure  is  the  rate  at  which  momentum  is  transferred  perpendicular  to  a  unit  surface  and 
for  an  ideal  gas  this  transfer  is  entirely  due  to  the  ballistic  motion  of  the  particles.  In  DSMC  there  is  also 
the  collisional  transfer  of  momentum  when  particles  within  a  cell  execute  a  collision.  As  discussed  in 
Section  2,  this  collisional  transfer  introduces  an  error  in  the  viscosity.  Interestingly,  the  pressure  is 
unaffected  because  of  a  symmetry  in  the  selection  of  collision  partners,  specifically,  particles  moving 
towards  each  other  are  as  likely  to  collide  as  those  moving  apart.  Since  the  relative  velocity  is  uncorrelated 
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with  the  relative  position  the  virial,  Avij  *  is  zero  on  average  and  thus  the  collisional  contribution  to  the 
pressure  is  zero. 


5.2  Consistent  Boltzmann  Algorithm 

The  first  DSMC  variant  for  a  non-ideal  gas  was  the  Consistent  Boltzmann  Algorithm  (CBA),  which 
models  a  hard-sphere  gas.  [21]  In  CBA  the  collision  step  is  modified  to  change  the  virial  in  order  to  yield 
the  contribution  to  the  equation  of  state  due  to  the  finite  volume  of  the  particles  The  CBA  algorithm  is 
identical  to  DSMC  with  a  minor  addition:  After  a  pair  of  particles  collides  the  particles  are  shifted  by  a 
displacement,  d ,  equal  in  magnitude  to  the  hard-sphere  diameter,  d ,  and  in  the  direction  of  the  change  of 
the  relative  velocity.  That  is, 


d 


~  tv  I 


This  displacement  (illustrated  in  Figure  5)  creates  a  non-zero  virial  that  correctly  captures  the  equation  of 
state  for  a  dense  hard-sphere  gas,  once  the  collision  rate  is  augmented  with  the  Enskog  Y-factor 
(enhancement  of  the  collision  rate  as  a  function  of  the  volume  fraction). 


After 


Figure  5:  Correlation  of  relative  position  and  relative  velocity  in  a  hard-sphere 
collision,  (left)  Displacement  of  particles  in  a  CBA  collision  (right) 


An  alternative  to  CBA  for  modelling  a  hard-sphere  gas  is  Enskog-DSMC  in  which  collision  partners  are 
restricted  to  being  a  diameter  apart  and  moving  towards  each  other.  [22,23]  While  this  modification  to  the 
collision  selection  lowers  the  computational  efficiency  it  has  the  advantage  of  having  a  theoretical 
foundation  based  on  the  Enskog  equation. 

5.3  Consistent  Universal  Boltzmann  Algorithm 

The  CBA  scheme  may  be  generalized  to  model  an  arbitrary  equation  of  state  by  making  the  collisional 
displacement  a  function  of  the  local  density  and  temperature.  [24]  In  essence,  we  “reverse  engineer”  the 
displacement  such  that  it  yields  the  desired  virial  and  thus  the  desired  equation  of  state.  This  variant,  called 
the  Consistent  Universal  Boltzmann  Algorithm  (CUBA),  has  been  demonstrated  with  the  van  der  Waals 
equation  of  state  and  shown  to  capture  the  correct  dynamics  even  in  the  condensed  phase.  [25]  Although  it 
is  surprising  that  a  simple  variant  of  DSMC  can  model  a  liquid,  as  the  density  increases  the  computational 
efficiency  decreases  because  of  the  increasing  collision  rate  so  CUBA  is  not  competitive  with  Molecular 
Dynamics  for  the  modelling  of  liquids.  Nevertheless,  for  moderately  dense  gases  it  remains  a  useful 
option. 
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5.4  Stochastic  Hard  Sphere  Dynamics 

As  discussed  in  the  earlier  sections,  DSMC  is  being  used  in  multi-scale  simulations  of  Brownian  motion 
since  it  is  an  efficient  scheme  for  modelling  the  solvent.  For  this  class  of  problems  it  is  often  desirable  to 
use  a  dense-gas  variant  in  order  to  reduce  the  compressibility.  However,  all  the  variants  discussed  above 
have  the  flaw  that  the  variance  of  density  fluctuations  are  inconsistent  with  the  compressibility  given  by 
the  equation  of  state.  For  example,  in  CBA  the  number  of  particles  in  a  volume  is  Poisson  distributed 
independent  of  the  density  so  the  density  fluctuations  have  the  same  variance  as  for  an  ideal  gas. 

Stochastic  Hard  Sphere  Dynamics  (SHSD),  another  dense-gas  variant  for  DSMC,  is  designed  to  produce 
the  correct  fluctuation  spectrum  consistent  with  the  equation  of  state.  [26]  The  collision  process  for  SHSD 
is  illustrated  in  Figure  6.  Particles  move  ballistically  between  collisions,  as  in  standard  DSMC.  When  two 
particles,  i  and  j,  are  less  than  a  diameter  apart  (| r  d  <  d)  there  is  a  probability  rate  (3%/<7)  vn  for  them  to 

collide,  where  vn  =  -Vy  *  fij/|fij|;  particles  only  collide  if  they  are  approaching  (i.e.,  if  vn  >  0).  If  the  collision 
is  accepted  then  it  is  evaluated  deterministically  as  if  the  particles  had  a  hard-sphere  diameter  of  ds  =  \r-f 
The  SHSD  model  has  dynamics  equivalent  to  a  fluid  with  a  linear  core  pair  potential.  [27] 


Figure  6:  Schematic  illustration  of  a  Stochastic  Hard  Sphere  Dynamics  (SHSD)  collision 


6.0  CONCLUDING  REMARKS 

Direct  Simulation  Monte  Carlo  is  for  the  computational  scientist  what  the  Boltzmann  equation  is  to  the 
mathematical  theorist.  The  popularity  of  DSMC  over  the  past  50  years  has  been,  in  large  part,  to  its 
intuitive  simplicity.  Readers  who  have  never  used  the  method  are  strongly  encouraged  to  set  aside  an 
afternoon  and  write  a  simple  DSMC  program  to  calculate  viscosity  (the  author’s  textbook  [3]  describes 
this  exercise  in  greater  detail).  Once  the  program  is  completed  (and  debugged!)  it  will  produce  the 
Chapman-Enskog  result  to  good  accuracy  in  about  ten  minutes.  When  you  contrast  this  with  the  tedious 
mathematical  calculation  required  using  the  Boltzmann  equation  you  will  appreciate  the  power  and  utility 
of  DSMC. 
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